perm filename SAIPOW.FAI[S,AIL]1 blob sn#102554 filedate 1974-05-22 generic text, type T, neo UTF8
COMPIL(POW,<FPOW,POW,LOGS,FLOGS,EXP$,LOG$>
	,<X11,X22,X33,OVPCWD>,<POW, FPOW, LOGS, FLOGS -- EXPON. ROUTINES>)
IFN ALWAYS,<	BEGIN	UTILS>
HERE(FPOW)
	SKIPA 	USER,-1(P)	;BASE
HERE(POW)
	FLOAT	USER,-1(P)
FPX:	MOVM	LPSA,-2(P)	;GET ABS(EXPONENT)
	JUMPE	LPSA,EXZERO	;0 EXPONENT
	MOVSI	A,(1.0)		;SET FOR FLOATING	
	JRST	2,@[FEXS]	;CLEAR AR FLAGS
FEXL:	ASH	LPSA,-1		;PREPARE TO LOOK AT NEXT BIT.
	FMPR	USER,USER	;SQUARE BASE
	 JFOV	 FPOWOV		;OVERFLOW/UNDERFLOW
FEXS:	TRZE 	LPSA,1		;COLLECT PRODUCT?
	FMPR	A,USER		;YES
	 JFOV	 FPOWOV		;OVERFLOW?
	JUMPN	LPSA,FEXL	;LOOP UNTIL EXPONENT ZERO.
	SKIPGE	-2(P)		;POSITIVE EXPONENT?
	   JRST	FEXDU1
POWRET: SUB	P,X33
	JRST 	@3(P)
FEXDU1:	MOVM	LPSA,A		;CHECK FOR OVERFLOW POSS.
	CAMGE	LPSA,[XWD 2400,1] ;SMALL NUMBER 
	 JRST	 FPDOV		;CALL UNDERFLOW
	MOVSI	LPSA,(1.0)	;TAKE RECIPROCAL OF ANS.
	FDVRM	LPSA,A
	JRST	POWRET		;AND RETURN IT.
EXZERO:	SKIPN	USER		;0↑0
ZRET:	 TDZA	 A,A		;RETURN 0
	MOVSI	A,(1.0)		;RETURN FLOATING 1
	JRST 	POWRET
FPOWOV:	SKIPN	TEMP,OVPCWD	;IF TRAPS ENABLED, USE EM
	 JSP	 TEMP,.+1	;ELSE GET FLAGS THIS WAY
	TLNE	TEMP,100	;SKIP IF NOT UNDERFLOW
FPDOV:	 MOVNS	 -2(P)		;UNDERFLOW -- CHANGE EXPONENT SIGN.
	MOVE	A,[XWD 400000,1] ;LARGE NEGATIVE NUMBER
	SKIPG	TEMP,-2(P)	;CHECK SIGN OF EXPONENT.
	 MOVEI	 A,0		;NEGATIVE ==> RESULT 0.
	SKIPGE	-1(P)		;CHECK SIGN OF BASE.
	 TRNN	 TEMP,1		;XOR SIGN OF EXPONENT.
	 MOVNS	 A		;MAKE +- LARGE NUMBER
	ERR	<Exponentiation under or overflow>,1
	JRST	POWRET		;RETURN.
HERE(FLOGS)
.FLOGS:	SKIPA	USER,-1(P)	;FLOATING BASE
HERE(LOGS)
.LOGS:	FLOAT	USER,-1(P)	;FLOAT THE BASE
	SKIPN	-2(P)		;IF ZERO EXPONENT,
	 JRST	 EXZERO		;GO TO COMMON CODE.
	MOVM	TEMP,-2(P)	;CHECK TO SEE IF 'FIX' WILL 
	CAMLE	TEMP,C1		;OVERFLOW
	 JRST	 USLGEP		;YES -- GO TO LOG-EXP
	FIX	TEMP,-2(P)	;CHECK TO SEE IF EXPONENT
	FLOAT	LPSA,TEMP	;HAPPENS TO BE AN INTEGER
	CAMN	LPSA,-2(P)	;IF SO, USE LOOPS TO
	 JRST	 [MOVEM TEMP,-2(P) ;BE SURE OF CORRECT SIGN
		  JRST FPX]
USLGEP:	JUMPE	USER,[SKIPGE -2(P) ;IF BASE ZERO, AND EXPT NEG.
			JRST FPDOV ;RETURN LARGE NUMBER
			JRST ZRET] ;ELSE RETURN ZERO.
	PUSH	P,USER		;ARGUMENT TO 'ALOG'
	PUSHJ	P,.LOG		;CALL IT.
	FMPR	A,-2(P)		;MULTIPLY BY EXPONENT
	PUSH	P,A		;ARGUMENT TO 'EXP'
	PUSHJ	P,.EXP		;CALCULATE
	JRST	POWRET		;AND RETURN.
C1:	243777777777		;2↑35 - EPSILON
HERE(EXP$)
.EXP:	PUSH	P,[0]		;ONE WORKING CELL
	PUSH	P,B		;AND ONE SAVED AC
	MOVE	LPSA,-3(P)	;GET ARGUMENT
	MOVM	A,LPSA		;GET ABSF(X)
	CAMG	A, E7		;IS ARGUMENT IN PROPER RANGE?
	JRST	EXP1		;YES, GO TO ALGORITHM
	ERR <EXP: under or overflow>,1
	HRLOI	A, 377777	;GET LARGEST FLOATING NUMBER
	SKIPG	LPSA		;WAS THE ARGUMENT POSITIVE?
	MOVEI	A, 0		;NO, RETURN 0
	JRST	EXPXIT		;AND RETURN
EXP1:	MULI	LPSA,400	;SEPARATE FRACTION AND EXPONENT
	TSC	LPSA,LPSA	;GET A POSITIVE EXPONENT
	MUL	TEMP,E5		;FIXED POINT MULTIPLY BY LOG2(E)
	ASHC	TEMP,-242(LPSA)	;SEPARATE FRACTION AND INTEGER
	AOSG	TEMP		;ALGORITHM CALLS FOR MULT. BY 2
	AOS	TEMP		;ADJUST IF FRACTION WAS NEGATIVE
	HRRM	TEMP,B 		;SAVE FOR FUTURE SCALING
	JUMPG	USER,ASHH	;GO AHEAD IF ARG GREATER THAN 0
	TRNN	USER,377	;ALL THESE BITS 0?
	 JRST	 ASHH		;YES -- GO AHEAD
	ADDI	USER,200	;NO -- FIX UP
ASHH:	ASH	USER, -10	;MAKE ROOM FOR EXPONENT
	TLC	USER, 200000	;PUT 200 IN EXPONENT BITS
	FADB	USER, -1(P) 	;NORMALIZE, RESULTS TO USER AND E
	FMP	USER,USER	;FORM X↑2
	MOVE	A, E2		;GET FIRST CONSTANT
	FMP	A, USER		;E2*X↑2 IN A
	FAD	USER, E4	;ADD E4 TO RESULTS IN USER
	MOVE	LPSA, E3	;PICK UP E3
	FDV	LPSA,USER	;CALCULATE E3/(F↑2 + E4)
	FSB	A,LPSA		;E2*F↑2-E3(F↑2 + E4)**-1
	MOVE	TEMP,-1(P)  	;GET F AGAIN
	FSB	A, TEMP		;SUBTRACT FROM PARTIAL SUM
	FAD	A, E1		;ADD IN E1
	FDVM	TEMP, A		;DIVIDE BY F
	FAD	A, E6		;ADD 0.5
	FSC	A, (B)		;SCALE THE RESULTS
EXPXIT:	POP	P,B		;RESTORE AC
	SUB	P,X33		;ADJUST STACK
	JRST	@2(P)		;RETURN.
E1:	204476430062		;9.95459578
E2:	174433723400		;0.03465735903
E3:	212464770715		;617.97226953
E4:	207535527022		;87.417497202
E5:	270524354513		;LOG(E), BASE 2
E6:	0.5
E7:	207540071260		;88.028
HERE(LOG$)
.LOG:
	SKIPGE	-1(P)		;CHECK SIGN OF ARGUMENT.
	ERR <LOG: Negative argument -- real part returned>,1
	MOVM	LPSA,-1(P)   	;GET ABSF(A)
	JUMPE	LPSA, LZERO	;CHECK FOR ZERO ARGUMENT
	CAMN	LPSA, ONE	;CHECK FOR 1.0 ARGUMENT
	JRST	ZERANS		;IT IS 1.0 RETURN ZERO ANS.
	ASHC	LPSA, -33	;SEPARATE FRACTION FROM EXPONENT
	ADDI	LPSA, 211000	;FLOAT THE EXPONENT AND MULT. BY 2
	MOVSM	LPSA,USER	;NUMBER NOW IN CORRECT FL. FORMAT
	MOVSI	LPSA, 567377	;SET UP -401.0 IN LPSA
	FADM	LPSA,USER 	;SUBTRACT 401 FROM EXP.*2
	ASH	TEMP, -10	;SHIFT FRACTION FOR FLOATING
	TLC	TEMP, 200000	;FLOAT THE FRACTION PART
	FAD	TEMP, L1	;TEMP = TEMP-SQRT(2.0)/2.0
	MOVE	LPSA,TEMP	;PUT RESULTS IN LPSA
	FAD	LPSA, L2	;LPSA = LPSA+SQRT(2.0)
	FDV	TEMP,LPSA	;TEMP = TEMP/LPSA
	MOVEM	TEMP,A		;STORE NEW VARIABLE IN A
	FMP	TEMP,TEMP	;CALCULATE Z↑2
	MOVE	LPSA, L3	;PICK UP FIRST CONSTANT
	FMP	LPSA,TEMP	;MULTIPLY BY Z↑2
	FAD	LPSA, L4	;ADD IN NEXT CONSTANT
	FMP	LPSA,TEMP	;MULTIPLY BY Z↑2
	FAD	LPSA, L5	;ADD IN NEXT CONSTANT
	FMP	A,LPSA		;MULTIPLY BY Z
	FAD	A,USER		;ADD IN EXPONENT TO FORM LOG2(X)
	FMP	A, L7		;MULTIPLY TO FORM LOGE(X)
LOGXIT:	SUB	P,X22
	JRST	@2(P)
LZERO:	ERR	<LOG: Argument 0; minus infinity returned>,1
	SKIPA	A, MIFI		;PICK UP MINUS INFINITY
ZERANS:	MOVEI	A,0		;MARG ANS ZERO
	JRST	LOGXIT		;AND RETURN
ONE:	201400000000
L1:	577225754146		;-0.707106781187
L2:	201552023632		;1.414213562374
L3:	200462532521		;0.5989786496
L4:	200754213604		;0.9614706323
L5:	202561251002		;2.8853912903
L7:	200542710300		;0.69314718056
MIFI:	400000000001		;LARGEST NEGATIVE FLOATING NUMBER
ENDCOM (POW)